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1. Introduction 

Light charged leptons constitute a particularly clean trigger for high-energy experiments. It 
is thus relatively easy to measure, despite their being small, the cross sections of processes 
that feature several final-state leptons. Prominent among such processes are those mediated 
by a pair of electroweak vector bosons which, depending on their identities, may give rise to 
a missing-energy signature as well. Vector boson pair production is interesting in at least 
two respects. Firstly, it is an irreducible background to Higgs signals, in particular through 
the W + W~ and ZZ channels which are relevant to searches for a standard model Higgs 
of mass larger than about 140 GeV. While always smaller than the W + W~ channel, ZZ 
decays may provide a cleaner signal due to the possibility of fully reconstructing the decay 
products of the two Z's. Secondly, di-boson cross sections are quite sensitive to violations 
of the gauge structure of the Standard Model, and hence are good probes of scenarios 
where new physics is heavy and not directly accessible at the LHC, yet the couplings in 
the vector boson sector are affected. 

The Higgs being such an elusive particle, searches require a combination of data-driven 
and theoretical methods in order to overcome irreducible backgrounds. From the theory 
viewpoint, this essentially implies computations as accurate as possible, and ideally able 
to realistically reproduce actual experimental events. On the other hand, modifications to 
tri-linear gauge couplings will tend to manifest themselves at large transverse momenta, 
and therefore this is the region where theoretical results have to be reliable. The inclusion 
of next-to-leading order (NLO) QCD corrections into four-lepton cross section predictions 
is a minimal, and fairly satisfactory, way to fulfill both of the two conditions above. 
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The aim of this paper is that of studying four-lepton hadroproduction at the NLO 
accuracy in QCD, by also adding the (finite) contribution due to gg fusion (which is for- 
mally of next-to-next-to- leading order, NNLO, yet enhanced by the gluon PDFs), and by 
including the matching with parton showers according to the MC@NLO formalism [1]; this 
is done by means of the framework aMC@NLO [2, 3] . Since aMCONLO provides the 
full automation of the matching procedure, and of the underlying matrix element compu- 
tations, we can give results for any four-lepton final state. In order to be definite, we have 
chosen to consider the process: 



with £ ,£' = e , \i. This choice is simply motivated by the fact that ZZ production as 
currently implemented in the MC@NLO package [4] does not feature production spin cor- 
relation and off-shell effects (which on the other hand are included, although in an ap- 
proximated way, in the case of W~W + and W^Z production in the codes of ref. [4]), and 
virtual-photon contributions with their interference with the Z's. The present aMC@NLO 
application remedies to these deficiencies by computing in an exact manner all relevant 
matrix elements (including the contributions of singly-resonant diagrams, which are po- 
tentially relevant to analyses in kinematical regions where one of the Z bosons is forced to 
be off-shell), and by including, as was already mentioned, the ^-initiated contribution as 
well which, although perturbatively suppressed, can be numerically important at the LHC 
owing to its dominant parton luminosity. We remind the reader that ZZ production at 
parton level and NLO accuracy in QCD has been studied for two decades now. The on-shell 
calculations of refs. [5, 6] have been subsequently improved to include leptonic decays [7], 
singly-resonant diagrams [8], and anomalous couplings [9]; these results have been used 
recently to match the NLO computation to parton showers according to the POWHEG 
formalism [10]. The 0{a 2 s ) process gg — > ZZ was first computed in refs. [11, 12, 13]. 
These papers have been later superseded by the inclusion of off-shell effects and Z/j* 
interference [14, 15]. As far as MC@NLO results are concerned, and on top of the phe- 
nomenological results presented here, this paper is a first for three reasons: 

• The matching with Pythia6 [16] for a kinematically non-trivial process. 

• The use of MadLoop [17] for the computation of loop-induced processes (i.e., the 
contributions of finite one- loop amplitudes squared). 

• The computations of scale and PDF uncertainties with a reweighting technique. 

In fact, the use of Pythia for the shower phase in MC@NLO has been limited so far to 
a proof-of-concept case [18], while the squares of loop amplitudes have been considered in 
MadLoop only for pointwise tests [17], and not in phenomenological applications. Finally, 
the reweighting procedure which we describe in this paper, and that is used to compute 
the scale and PDF dependences of the cross section, allows one to determine the corre- 
sponding uncertainties without requiring any additional computing time. We point out 
that the capabilities listed in the three items above are not process dependent and are fully 
automated. 




(1.1) 
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This paper is organized as follows: in sect. 2 we discuss some of the technicalities 
relevant to four-lepton production, and the procedure implemented in aMC@NLO for the 
determination of scale and PDF uncertainties. In sect. 3 we present selected phenomeno- 
logical results, and in sect. 4 we draw our conclusions. Further details on scale and PDF 
dependence computations are reported in appendix A. 

2. Technical aspects of the computation 

The framework of the aMC@NLO programme used for the phenomenology studies of 
this paper is unchanged w.r.t. that employed in ref. [3]; in particular, the underlying 
tree-level computations are performed with MadGraph v4 [19]. We remind the reader 
that aMC@NLO automates all aspects of an NLO computation and of its matching with 
partons showers. One-loop amplitudes are evaluated with MadLoop [17], whose core 
is the OPP integrand reduction method [20] as implemented in CutTools [21]. Real 
contributions and the corresponding phase-space subtractions, achieved by means of the 
FKS formalism [22], as well as their combination with the one-loop and Born results and 
their subsequent integration, are performed by MadFKS [23], which finally takes care of 
the MC@NLO matching [1] as well. 

The novel features of aMC@NLO whose results we present here are the matching with 
the virtuality-ordered Pythia6 shower, based on ref. [18]; the possibility of determining 
through reweighting the scale and PDF uncertainties affecting our predictions, which we 
shall discuss in general in sect. 2.1 (with some further technical details given in appendix A); 
and the use of MadLoop for the computation of the 0(a 2 s ) partonic process 



The matrix elements relevant to this process are UV- and IR-finite, and therefore the 
corresponding generated events are treated as an LO sample from the viewpoint of matching 
with parton showers. The amplitudes can be straightforwardly computed by MadLoop. In 
the current version, MadLoop assumes that these amplitudes will have to be multiplied by 
Born ones; we have extented the scope of the code, in order for it to compute amplitudes 
squared. 

The fact that the process in eq. (2.1) is not a genuine virtual correction, i.e. it lacks an 
underlying Born, implies that event unweighting is essentially a brute-force operation (as 
opposed to the normal situation where the presence of an underlying Born allows one to 
pre-determine with good accuracy the peak structure of the one-loop matrix elements, thus 
significantly increasing the unweighting efficiency). This results in a fairly large number of 
calls to the matrix elements of eq. (2.1) per unweighted event and, given the computing 
performances of the current version of MadLoop, renders it very time expensive to obtain 
a good-sized unweighted-event sample (say, 0(1 M)). We have therefore opted for an 
alternative approach. We have generated unweighted events using the matrix elements of 
the Born process 




(2.1) 
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multiplied by the gg parton luminosity (rather than by the uu one that would be required if 
computing the Born contribution proper). The kinematical configurations so obtained have 
been used to compute the matrix elements of eq. (2.1), which thus provide the correct event 
weights. A similar procedure has been employed in a recent study on Higgs production via 
heavy-quark loops and has been shown to work well even for more complex final states [24] . 

2.1 Scale and PDF uncertainties 

Among the parameters which enter a short-distance cross section, renormalization ((J, R ) and 
factorization (fi F ) scales and PDFs play a special role, since their variations are typically 
associated with the purely theoretical uncertainty affecting observable predictions (which, 
in the case of PDFs, is not quite correct - one may say that PDF variations parametrize the 
uncertainties not directly arising from the process under study). It is therefore fortunate 
that (the bulk of) the time-consuming matrix element computations can be rendered inde- 
pendent of scales and PDFs, as opposed to what happens in the case of other parameters, 
e.g. particle masses. This is doable thanks to the fact that short-distance cross sections 
can be written as linear combinations of scale- and PDF-dependent terms, with coefficients 
independent of both scales and PDFs; it is thus possible to compute such coefficients once 
and for all, and to combine them at a later stage with different scales and PDFs at essen- 
tially zero cost from the CPU viewpoint. The crucial (and non-trivial) point is that the 
noteworthy structure mentioned before is a feature not only of the parton-level LO and 
NLO cross sections, but also of the MC@NLO ones. This implies that from the conceptual 
point of view the same procedure for determining scale and PDF uncertaintities can be 
adopted in MC@NLO as in LO-based Monte Carlo simulations; for the former we shall 
simply need to compute a larger number of coefficients than for the latter. 

In order to illustrate what is done in aMC@NLO, we start from dealing with LO 
computations, for which the notation is simpler. Here and in what follows, the expressions 
for the short-distance cross sections are taken from ref. [23]. The fully-differential cross 
section is: 

da (LO) = /i(xi, fi F )f 2 (x2, ii F )da^J B jd X Bj , (2.3) 

d(J (B,n) = M (n,0)±_ d(j)n _ (2 _ 4) 

Here, d\Bj denotes the integration measure over the Bjorken x's, and Jsj is the (possibly 
trivial) corresponding jacobian factor. The quantities on the r.h.s. of eq. (2.4) are the 
scattering amplitude squared J\A^ n ' ^ (including coupling constants, and colour/spin average 

(S) 

and flux factors), a set of cuts J n L which prevent phase-space singularities from appearing 
at this perturbative order, an average factor ftf for identical final-state particles, and the 
n-body phase space d(f> n . We shall write the latter as follows: 

#„ = * (n) (MXn))dXn, (2-5) 

where d\ n is a measure that understands the choice of 3n — 4 independent integration 
variables, and fC n (x n ) denotes the four- momentum configuration associated with a given 
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choice of the latter. Finally, <&( n ^ is the phase-space factor (including jacobian) that arises 
from the choice of d\n- One can then rewrite eq. (2.3) as follows: 

da^ = h{x l ,ii F )h{x 2 ,ii F )gf( j i R )w^d X B 3 dxn, (2.6) 

with 

w (B,n) _ M<- n ® J n ^ ) (n) 

and, consistently with ref. [23], b is the power of a s implicit in Ai^ n '°\ By construction, 
the quantity w^ B ' n ^ defined in eq. (2.7) is independent of scales and PDFs, and is the first 
example of the coefficients mentioned at the beginning of this section. 
When integrating eq. (2.3), one obtains the iV-event set 

|/C„ ; j , x\-i , X2-i , Sj | , (2.8) 

with 

da (LO) 

Ei = — ()C n .i, xi-i, x 2; i) , = constant , (2.9) 

dXBjdXn 

in the case of weighted and unweighted event generation respectively 1 . Given an observable 
O, the cross section in the range O low < O < VPP (e.g., an histogram bin) will be 

N 

e(o (/c n;i ) - o LOW ) e(o UPP - o (/c n;i )) (2.10) 

i=i 

at the hard-subprocess level, and 

N 

©(O (MC(/C n;i , x 2]l )) - O low ) e(0 UPP - O (MC(/C n;i , x 2;i ))) Hj (2.11) 

i=i 

after shower. In eq. (2.11), MC()C n -i,xi-i,x 2 -i) denotes the complete final-state configura- 
tion obtained by showering the hard subprocess. 

When making alternative choices for scales and PDFs, say fi' R , fi' F , and f-, one can 
simply repeat the procedure outlined above; this is straightforward, but extremely time 
consuming if the full scale and PDF uncertainties have to be determined. Alternatively, 
one can compute eqs. (2.10) and (2.11) by performing the replacement 

El — > EiTZi, (2.12) 

where 

/<i£r (LO) 

-, -7- (]Cn;i,X 1 . i ,X 2 . i ). (2.13) 
dXBjdXn 



x It should be obvious that the kinematic configurations and Bjorken x's in eq. (2.8) are not the same 
for weighted- and unweighted-event generation. Nevertheless, we have used an identical notation since no 
confusion is possible. 
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The denominator here is computed using eq. (2.6), i.e. with the original choices of scales 
and PDFs. As we shall show in the following, an analogous equation will hold for both the 
NLO and MC@NLO cases. From eq. (2.6) we obtain: 



fl(xi;i, ^F)f2(x2;i,^F)gf(HR) 



which shows explicitly that the computation of IZi does not entail any matrix-element 
calculation. Therefore, after performing the bulk of the calculation, i.e. the determination 
of the set in eq. (2.8), one can fill a "central" histogram with the weights defined in 
eq. (2.9), and as many "variation" histograms as one likes with the weights that appear 
on the r.h.s. of eq. (2.12), by simply recomputing the factors IZi for all the different scales 
and PDFs needed. 

The procedure outlined above is exact when applied to eq. (2.10), but it implies an 
approximation in the case of eq. (2.11). This is because the quantity MC(/C n; j, xi-i, X2-i) 
does contain an implicit dependence on /j , and these functions are not replaced by j[ when 
using eq. (2.12) (since such a replacement would imply performing the shower for each 
new choice of PDFs). However, this approximation is usually a very good one (barring 
perhaps the corners of the phase space, and in particular the large-rapidity regions). This 
is empirically well known, and is due to the fact that the Sudakov form factors used in 
the backward evolution of initial-state partons are sensitive to the ratios of PDFs, whose 
variation is much smaller than that of their absolute values (the more so within a PDF- 
error set, which is the typical application of the procedure discussed here). Furthermore, it 
should be clear that both scale and PDF uncertainties are defined using well-motivated but 
ultimately arbitrary conventions, and therefore are not quantities that can be computed 
with arbitrary precision. 

We now turn to discussing the case of NLO short-distance cross sections. Equa- 
tion (2.6) is generalized as follows: 

d(J (NLO) = ^ d<7 (NLO ' a) , (2.15) 

a 

da ^o, a) = /l(x W,^) )/2 (4"),^))^W dxB ^ Xn+1 , (2.16) 

where a = E, S, C, and SC correspond to the contributions of the fully-resolved config- 
uration (the event), and of its soft, collinear, and soft-collinear limits (the counterevents) 
respectively. The quantities will be written as follows: 



= gf + \^ 



( (")\ 2 / 



Q J \ Q J 

+ 9f^n')W B S aS , (2.17) 

where Q is the Ellis-Sexton scale; in eq. (2.17), this scale (which is extensively used in the 
manipulation of the one-loop contribution - see ref. [23] for a discussion in the context 
of MadFKS) offers a convenient way to parametrize the dependence on the factorization 
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and renormalization scales. As the notation in eqs. (2.16) and (2.17) suggests, the values 
that such scales assume in the event and counterevents may be different from each other 
(however, they must tend to the same value when considering the relevant infrared limits). 
Furthermore, eq. (2.16) takes into account the fact that, when using event projection 
to write an NLO partonic cross section (see ref. [1]), the Bjorken x's of the event and 
counterevents need not coincide. The last term on the r.h.s. of eq. (2.17) is the Born 
contribution, which in MadFKS can be integrated simultaneously with the other terms 
that enter an NLO cross section; having a soft-type kinematics, it is naturally associated 
with the soft counterevent. The coefficients Wq U \ W F a \ and introduced in eq. (2.17) 
are the analogues of u>( B > ra ) defined in eq. (2.7) - they are scale- and PDF- independent; 
their explicit forms are given in appendix A. As far as W B is concerned, it is equal to 
W ( B ' n \ up to a normalization factor (due to the different integration measures in eq. (2.6) 
and (2.16)). 

The integration of the NLO cross section leads to the set of N weighted events 2 : 




«) J a ) J a ) L I ( 2 IS) 

> X l;i > x 2;i ) "i L-rc^e^l ' 1 ' 

1=1 



with 

Equation (2.10) is generalized as follows: 

N 

EE e (° - K p - ° ) - a) • (2 - 20) 

i=l a 

As is well known, for any given % the weights defined in eq. (2.19) may diverge, but the sum 
in eq. (2.20) is finite for any infrared-safe observable. When changing scales and PDFs, we 
can now adopt the same procedure introduced before for the LO cross sections, and define 
the rescaling factors: 

( Q ) _ fi(J a ) ,A a )\f>(„( a ) ,/(«hJ 



' C i — Jl\ X l-A i^F )j2{ X 2;i'f i 



Wt ] (4 Q i,) + W F a) (JCi% t ) log (^j + W**> (4 Q i,) log 



+ 9f^V)W B {K^,.)5 a ^ I ^(Xgy,^,^) • (2.21) 

The quantities 7?.j a ' > thus defined do not factorize the coefficients W , and hence have a 
more complicated form than their LO counterpart, eq. (2.14). Nevertheless, precisely as in 



2 We remind the reader that parton-level, NLO events not matched with showers cannot be unweighted 
without the introduction of an arbitrary and unphysical cutoff. 
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the case of eq. (2.14), their computations only require the evaluation of PDFs and coupling 
constants for each new choice of PDFs and scales, the idea being that the coefficients W 
are computed once and for all when integrating the partonic cross section, and are stored 
for their later re-use in eq. (2.21). 

We have also understood that the formulae presented above apply to the contribution 
due to one given FKS pair, i.e. they are proportional to a given S function Sij. We remind 
the reader that while the observable cross section is obtained by summing over all S- 
function contributions, these are fully independent from each other, and can be effectively 
treated as independent partonic processes. 

We finally discuss the case of MC@NLO. The short-distance cross sections can be 
written as follows [1]: 

da m = da (Nh0 ' E> - £ da (MC ' c) , (2.22) 

c 

da (S) = ^Ar (MC ' c) + da (Kh °> a \ (2.23) 

c a=S,C,SC 
with the Monte Carlo counterterms given by: 

da^ = h(x^ c \ Mf )/ 2 (4 mc ' c) , fi F )gf +2 (ti R )w^d XBj d Xn+1 . (2.24) 

The index c in eqs. (2.22) and (2.23) runs over all particles which are colour-connected with 
the particle that branches. The latter is identified with the parent particle of the FKS pair 
which determines the S function Sij implicit in <icT (NLO ' Q) . This is sufficient to render 
eqs. (2.22) and (2.23) locally finite, as prescribed by the MC@NLO formalism. Note that 
the sum over FKS pairs is equivalent to summing over all possible MC branchings. By 
integrating eqs. (2.22) and (2.23) one obtains the sets of events: 

whose kinematic parts are: 

£m = {/c H;i {*S C < C \ *g c < c) } c } , (2.26) 

£s-,i = {JC s ,Ax$,^}«=s,c,scAx^ C \x^ C) }^ ■ (2-27) 
The final-state kinematic configurations relevant to H and § events are: 

^e;j = ^n+l-i , K-S;i = > (2.28) 

where we have used the fact that in MadFKS and aMC@NLO the phase-space parametriza- 
tion is chosen in such a way that the kinematics of the counterevents coincide (up to the 
irrelevant unresolved partons; see ref. [23]). The weights that appear in eq. (2.25) are: 

aXBjdXn+i aXBjdXn+i 
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for weighted events, and 



n(H) 




n(S) 









constant 



(2.30) 



for unweighted events. One can now define rescaling factors analogous to those of eqs. (2.13) 
and (2.21): 



(H) 



W : o (£) fe) + wf ) (f H;i )logi 



+ W ? f ) fe)log 



'AT 



- E /i(*S°' c) > 4 E) )/^(4T' C) ' /4* V + Vi E V (MC ' c) 



dXBjdXn 



+1 



■(^H;j) , 



(2.31) 



ft. 



(S) 



a>=S,C,SC 



a f +2 (A s) ) | whsm) + ^r(%) log ^ 
v., ,as) - ... . i / ^ (s) 



/4f 

Q , 



, Q . 



+ 9f^T')W B (£s;i)S 



(2.32) 



/ dxBjdXn+i 

Several observations are in order here. Firstly, the weights defined in eq. (2.29) are finite, 
at variance with those relevant to the computation of the parton-level NLO cross section, 
eq. (2.19). This is due to the properties of the MC counterterms, and to the fact that 
the kinematics of the H and S events is uniquely defined. Secondly, and as a consequence 
of the former point, the quantities relevant to the counterevents are summed together in 
eq. (2.32), at variance with what happens in eq. (2.21). Thirdly, in eqs. (2.31) and (2.32) 
the scales entering the counterevents have been set equal to the values assumed in the 
soft configurations, and those entering the MC counterterms equal to the corresponding 
NLO event contributions. Many variants of these choices are possible (which is part of the 
ambiguity in the definition of the scale dependence in an NLO computation), and those 
given here are simply the current defaults in aMC@NLO. 

Using the definitions given above, eq. (2.11) is generalized as follows: 



£ 

i=i 



0(0 (MC(£ Hi i)) - Olow) ©(Oupp - O (MC(5 H;j ))) sf } nf> 
E ©(O (MC(£ S ;i)) - Olow) ©(O upp - O (MC^))) nf ] . 



(2.33) 
(2.34) 
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/i(si)Born , (2.35) 



Here, the same remark as made in the case of LO-showered predictions applies, that con- 
cerns the PDF dependence implicit in MC(%) and MC(£§ ; j). There is however a further 
subtlety in the case of MC@NLO. When showering the Born contribution to § events, one 
cancels the MC counterterm contribution to H events. In the case of initial-state emis- 
sions, we can schematically write the contribution to the H events obtained by showering 
the Born as follows: 

h{xi/z) 
fi(xi) 

where we have assumed an emission from leg 1 to be definite, and all overall factors ir- 
relevant to this discussion have been neglected. The ratio of PDFs in eq. (2.35) is due to 
the MC branching, and [/i(xi)Born] is the short-distance contribution. When S events are 
reweighted with 1Z\ without changing the PDFs in the shower, eq. (2.35) becomes: 

ta^Born] , (2.36) 
Jl{ x l) L J 

and therefore the cancellation of the PDFs with argument x\ does not take place any longer. 
Since changing the PDFs in the shower would amount to replacing f\ with /{ in eq. (2.36), 
eq. (2.36) implies a fractional difference w.r.t. the use of /{ in both the short-distance and 
shower computations equal to: 

h{xi/z) f[{xi) 

f>( Xl /z) ' {Z - 6 '> 

which is a number generally fairly close to one. In fact, what is discussed here is precisely 
the situation one faces when running MC@NLO and adopting different PDFs in the short- 
distance calculation and shower phases. This is rather commonly done, and no significant 
differences are observed w.r.t. the results obtained with the strictly correct procedure of 
using the same PDFs in the two phases. While this approximation may be undesirable 
when very accurate results are mandatory, it is always acceptable when computing the 
PDF uncertainty, which is by definition not a precision result. Furthermore, it should be 
kept in mind that the procedure currently used for determining PDF "error" sets is that of 
employing parton-level NLO cross sections; hence, such PDFs are not guaranteed to serve 
their purpose for those observables for which matched predictions differ significantly from 
the NLO ones (i.e., where resummation effects are important). 



3. Phenomenological results 

In this section we present some sample results obtained with aMC@NLO, its LO counter- 
part (aMC@LO), and for parton-level NLO (0(a s )) and LO (O(a s )) cross sections. We 
also give showered predictions for the gg 0{a 2 s ) contribution. In both HERWIG [25] and 
Pythia, we have switched QED showers off. We consider the production of the following 
leptonic final states: 

£+£-£'+£'- = e+e-fi+n' , (3.1) 
£+£-£+£- = e+e-e+e" , (3.2) 
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Parameter 


value 


Parameter 


value 


m z 


91.118 


r z 


2.4414 


a- 1 


132.50698 




1.16639-10" 5 


(LO), \ 

«s {mz) 


0.13939 


(NLO) , \ 


0.12018 



Table 1: Settings of physical parameters used in this work, with dimensionful quantities given in 
GeV. 

at the 7 TeV LHC. The only cut at the generation level is: 

M(^W^) > 30 GeV , (3.3) 

which is applied to all unlike-sign lepton pairs. This cut, that serves the purpose of avoiding 
the zero-virtuality divergences of diagrams with resonant photons, is obviously not neces- 
sary in the case of eq. (3.1) for pairs, but it is nonetheless applied there in order to 
perform a comparison on equal footing between the e + e~ [i + and e + e~e + e~ final states. 
The parameters used in this paper are reported in table 1, where the values of a s have 
been chosen as prescribed by the PDF sets we have adopted, MSTW2008(n)lo68cl [26] 
(with their associated "error" sets). The default renormalization and factorization scales 
are set equal to the invariant mass of the four leptons, 

Mo = M{t + rt^ + t^~) . (3.4) 

When studying scale dependences we shall vary fj, R and [i F independently, i.e. we shall 
set (fi R ,fj, F ) = (k r ho, k f [j>o) and consider the combinations (k r ,k f ) = (1,1), (1/2,1/2), 
(1/2, 1), (1, 1/2), (1,2), (2, 1), (2,2). The overall scale uncertainty will be the envelope of 
all individual results. 

We begin by presenting in table 2 the fully inclusive cross sections (bar the invariant- 
mass cuts of eq. (3.3)). For the NLO results and the gg contribution we also give the 
scale (first error) and PDF (second error) uncertainties. The scale dependence at the 
NLO is fairly small, which is in part due to the fact that the qq parton luminosity is 
almost /^-independent in the Bjorken-x range that gives the dominant contribution to 
the cross section. This implies that at the LO (which, being 0(a®), is independent of 
H R ) the scale uncertainty is comparable or smaller, and thus not particularly meaningful: 
we refrain from quoting it explicitly here. The PDF uncertainty is estimated following 
the prescription given by the PDF set [26] (asymmetric Hessian); at the LO, we obtain 
a fractional uncertainty similar to that of the NLO result. On the other hand, the gg 
channel displays the typical scale dependence of an 0{a 2 s ) LO cross section, the figures 
in the table receiving the dominant contributions from the /j, r dependence. The scale 
uncertainty is largely dominant over the PDF one. Not surprisingly, the ratio between the 
two leptonic channels considered in this paper is equal to two up to a very small correction 
(only slightly larger than the integration uncertainty), since differences arise in the off- 
resonance regions (where singly-resonant diagrams are more important than elsewhere) 
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Cross section (fb) 



Process 



qq/qg channels 



gg channel 



O(a s ) + 0(a s ) 




pp — > e^e n^fi 



9.19 



12.90 



+0.27(2.1%)+0.26(2.0%) 
-0.23(1.8%)-0.22(1.7%) 



0.566 



+0.162(28.6%)+0.012(2.1%) 
'-0.118(20.8%)-0.014(2.5%) 



pp — > e^e e^e 



4.58 



6.43 



+0.13(2. 1%)+0.11(1. 7%) 
'-0.13(2.0%)-0.10(1.6%) 



Table 2: Total cross sections for e + e~^ + /i~ and e + e~e + e~ production at the LHC (y/S = 7 TcV), 
for the production channels considered in this paper, and within the cuts given in eq. (3.3). The first 
and second errors affecting the results are the scale and PDF uncertainties (also given as fractions 
of the central values). See the text for details. 

whose contribution to the total integral is small. Given that the gg matrix elements do 
feature only doubly-resonant diagrams, we have not computed the result relevant to the 
e + e~e + e~ channel. Finally, we point out that the entries of table 2 have been compared 
with their analogues obtained from MCFM (also for the non-central scale values, which is 
a check of the reweighting technique discussed in sect. 2.1), and complete agreement has 
been found. 

As was mentioned above, the scale uncertainties reported in table 2 are the envelopes 
of the results obtained with the seven combinations of (fi R ,(i F ) considered here. It is 
obvious that this is not equivalent to variations in the range /xo/2 < n R ,fi F < if the 
dependence on either scale is not monotonic. While the study of continuous variations 
is clearly impossible, by means of reweighting one can probe the relevant scale ranges 
with more than three values for each scale, at a very modest CPU cost. In fact, each 
choice of (fi R ,fi F ) corresponds to filling an extra set of histograms (each set containing all 
observables of interest), with the same kinematics as for the central set and with rescaled 
weights, which is a negligible fraction of the generation-plus-shower computing time. 

We now turn to discussing the case of differential distributions, and we start by consid- 
ering the 0(pP s ) and 0(a s ) contributions. We shall later address the (generally very small) 
differences between HERWIG and Pythia6, and here we shall use only HERWIG in or- 
der to be definite. The aim of figs. 1-3 is that of assessing the perturbative behaviour of our 
predictions. In the main frame, these figures present the results for the e + e~[i + [i~ channel 
obtained with aMC@NLO (solid black), NLO (dashed red), aMC@LO (solid blue), LO 
(dashed magenta). The middle insets show the scale (dashed red) and PDF (solid black) 
uncertainties, both with aMC@NLO, and given as ratios over the central values. Finally, 
the lower insets present the ratios 



with O any of the observables considered, as computed by aMC@NLO. 

Figure 1 displays two observables constructed with the sum of the four-momenta of the 
four leptons, the invariant mass (left panel) and the transverse momentum (right panel). 




(3.5) 
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200 400 600 800 1 2 

M(e + e"M"V) [GeV] log 10 (p T (e + e>V")/GeV) 



Figure 1: Four-lepton invariant mass (left panel) and transverse momentum (right panel), as pre- 
dicted by aMC@NLO (solid black), aMC@LO (solid blue), and at the (parton-level) NLO (dashed 
red) and LO (dashed magenta). The middle insets show the aMC@NLO scale (dashed red) and 
PDF (black solid) fractional uncertainties, and the lower insets the ratio of the two leptonic channels, 
eq. (3.5). See the text for details. 

These have very different behaviours w.r.t. the extra radiation provided by the parton 
shower, with the former being (almost) completely insensitive to it, and the latter (almost) 
maximally sensitive to it. In fact, the predictions for the invariant mass are basically 
independent of the shower, with NLO (LO) being equal to aMC@NLO (aMC@LO) over 
the whole range considered. The NLO corrections amount largely to an overall rescaling, 
with a very minimal tendency to harden the spectrum. The four-lepton p T , on the other 
hand, is a well known example of an observable whose distribution at the parton-level LO 
is a delta function (in this case, at p T = 0). Radiation, be it through either showering or 
hard emission provided by real matrix elements in the NLO computation, fills the phase 
space with radically different characteristics, aMC@LO being meaningful at small p T and 
NLO parton level at large p T - aMC@NLO correctly interpolates between the two. The 
different behaviours under extra radiation of the two observables shown in fig. 1 is reflected 
in the scale uncertainty: while in the case of the invariant mass the band becomes very 
marginally wider towards large M(e + e~ p, + fj,~) values, the corresponding effect is dramatic 
in the case of the transverse momentum. This is easy to understand from the purely 
perturbative point of view, and is due to the fact that, in spite of being 0(a s ) for any 
p T > 0, the transverse momentum in this range is effectively an LO observable (the NLO 
effects being confined to p T = 0). The matching with shower blurs this picture, and in 
particular it gives rise to the counterintuitive result where the scale dependence increases, 
rather than decreasing, when moving towards large p T [18]. Finally, the lower insets of 
fig. 1 display the ratio defined in eq. (3.5) which, in agreement with the results of table 2, 
is equal to one half in the whole kinematic ranges considered. The only exception is the 
small invariant mass region, where off-resonance effects become relevant. 
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Figure 2: Same as in fig. 1, for the inclusive p T of the positively-charged leptons (left panel), and 
the inclusive p T of the same-charge lepton pairs (right panel), both with Z-id cuts. 



A gauge-invariant way to suppress off-resonance effects, and to select doubly-resonant 
contributions, is that of imposing: 

\M (£+£-) - m z \ < 10 GeV (3.6) 

on all equal-flavour lepton pairs; we call the cut of eq. (3.6) the Z-id cut. Lepton pairs that 
pass the Z-id cut are called Z-id matched, and can be roughly seen as coming from the 
decay of a (generally off-shell) Z boson. While in the case of the e + e~ fi + fi~ channel there 
is only one way to choose two same- flavour lepton pairs, there are two different pairings in 
e + e~e + e~~ production. In the case both of these pairings result in lepton pairs that fulfill 
eq. (3.6), we choose that with the smallest pair invariant mass, and assign the Z-id matched 
pairs according to this choice; in practice, this is a rare event. By imposing the Z-id cuts 
the M(£ + £~£^ + £^~) distribution falls steeply below threshold and gets no contributions 
below 160 GeV. 

In fig. 2 we present two transverse momentum distributions, relevant to the positively- 
charged leptons (left panel), and to same-charge lepton pairs (right panel); hence, there are 
two entries in each histogram for any given event. These results are obtained by applying 
the Z-id cuts, but we have in fact verified that without such cuts we obtain exactly the 
same patterns. In the case of the p T of the individual lepton, the aMC@NLO (aMC@LO) 
prediction is fairly close to the NLO (LO) one, but tends to be slightly harder, owing to 
the extra radiation generated by the shower. This effect is more pronounced at the LO 
than at the NLO, which is the sign of a behaviour consistent with perturbation theory 
expectations. In fact, at the LO all hadronic transverse momentum is provided by the 
shower, while at the NLO this is not the case; therefore, at the NLO the shower will have 
less necessity to "correct" the prediction obtained at the parton level, a tendency which is 
naturally embedded in a matching prescription such as aMC@NLO. The scale dependence 
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Figure 3: As in fig. 1, for the inclusive 77 of the opposite-charge, Z-id matched lepton pairs (left 
panel), and the inclusive distance of the opposite-charge, non-Z-id matched lepton pairs (right 
panel). 

is quite small over the whole range in p T , but tends to grow larger towards larger p T . This 
effect has the same origin as that observed in the right panel of fig. 1, but it is much more 
moderate than there. This is due to the fact that in the present case the whole range in p T 
is associated with complete NLO corrections. The PDF uncertainty is seen to be similar to 
or slightly smaller than that due to scale variation; parton densities are well determined in 
the x range probed here. Finally, there is no difference between the two leptonic channels 
for this observable; as already mentioned above, this conclusion is independent of whether 
one applies the Z-id cuts. The p T of the lepton pairs shown in the right panel of fig. 2 
follows the same pattern as the one we have just discussed, but the differences between 
the various predictions are larger in this case. In particular, aMC@LO is closer to NLO 
than to LO, which is a consequence of the more important role played by extra radiation in 
this case (as one expects, the present one being a correlation between two particles rather 
than a single-inclusive observable). Again, the closeness of NLO and aMC@NLO results 
shows the desired perturbative behaviour. The more significant impact of extra radiation 
on this variable is reflected in the slightly larger scale dependence at large p T 's w.r.t. what 
happens for the transverse momentum of the individual leptons discussed before. The two 
leptonic channels agree well, also when removing the Z-id cuts. 

Figure 3 shows two observables constructed after applying the Z-id cuts, namely the 
pseudorapidity of lepton pairs with opposite charge which are also Z-id matched (left 
panel; this is then the pseudorapidity of would-be Z bosons), and the azimuthal distance 
between leptons of opposite charge which are not Z-id matched (right panel; thus, these 
are leptons emerging from different would-be Z bosons). As in the case of fig. 2, there are 
two entries in each histogram for any given event. These two observables are dominated 
by small transverse momenta, and therefore it is not suprising that, at both 0(a®) and 
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M(e + e"M"V) [GeV] log 10 (p T (e + e>V")/GeV) 



Figure 4: Same observables as in fig. 1, for aMC@NLO+gg HERWIG (solid black) and Pythia 
(dashed blue) results. The rescaled gg contributions with HERWIG (open black boxes) and Pythia 
(open blue circles) are shown separately. Middle insets: scale (dashed red) and PDF (solid black) 
fractional uncertainties. Lower insets: aMC@LO/(aMC@NLO+.g3) with HERWIG (solid black) 
and Pythia (dashed blue). 



0(a s ), the predictions are quite independent of whether a shower is generated or not. 
Slight differences can be seen in the case of the A(p distribution, which is indeed known to 
be more sensitive than pseudorapidity to extra radiation. The small-p T dominance ensures 
that scale and PDF uncertainties are flat over the whole kinematic ranges, and of the order 
of those relevant to total cross section. 

We now discuss the impact of the 0{a 2 s ) gg channel on our predictions. The argument 
for considering such a channel, despite its being of the same perturbative order as all other 
NNLO contributions which cannot be included, is the dominance of its parton luminosity 
over those of the qq and qg channels. This dominance grows stronger with decreasing 
final-state invariant masses, and hence the 0(a 2 s ) versus NLO comparison is significantly 
influenced by the cut in eq. (3.3) - by lowering such a cut, the relative importance of the 
gg contribution will grow bigger than the 5%-ish reported in table 2. We also discuss in the 
following the differences that arise when matching our calculation to Pythia6 rather than 
to HERWIG. We remind the reader that, depending on input parameters, Pythia is rather 
effective in producing radiation in the whole kinematically-accessible phase space. This is 
not particularly useful in the context of a matched computation, where hard radiation 
is provided (in a way fully consistent with perturbation theory) by the underlying real- 
emission matrix elements. Therefore, we have set the maximum virtuality in Pythia 
equal to the four-lepton invariant mass. For consistency, this setting has been used also 
when showering the gg-initiated contribution. 

Figures 4, 5 and 6 present the same observables as figs. 1, 2 and 3 respectively. In 
the main frame, we show the aMC@NLO predictions plus the gg contribution (including 
shower), as resulting from HERWIG (solid black) and Pythia (dashed blue) - we shall 
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Figure 5: Same as in fig. 4, for the observables introduced in fig. 2. 



call these predictions aMC@NLO+g<? for brevity. The dashed-plus-symbol histograms are 
the gg contributions (open black boxes: HERWIG; open blue circles: Pythia), rescaled 
in order to fit on the same scale as aMC@NLO+<7g; the rescaling factors are equal either 
to 10 or 20, as indicated in the figures - they are chosen so as to render the comparisons of 
shapes as transparent as possible. The middle insets show the fractional scale (dashed red) 
and PDF (black solid) uncertainties of the gg contribution alone (i.e., the extrema divided 
by the central gg prediction); those shown are computed with HERWIG, but Pythia 
gives identical results. The lower insets present the ratios aMC@NLO/(aMC@NLO+gg) 
for HERWIG (solid black) and Pythia (dashed blue). 

The HERWIG and Pythia results for the four-lepton invariant mass shown in the 
left panel of fig. 4 are identical in the case of the gg contribution, and practically so in 
the case of aMC@NLO. This is because both event generators fix the invariant mass of 
the final-state system when doing initial-state showers; hence, any difference in the four- 
lepton invariant mass can only arise from the real-emission contribution to aMC@NLO, 
where the final-state system does not coincide with the four leptons. The gg channel 
has a softer distribution than the qq + qg ones, since by moving towards large masses one 
probes larger values of the Bjorken x's, where quark parton luminosities are relatively more 
important than at threshold. In the small-invariant-mass region, the gg contribution does 
not exhibit any peak structure, owing to the lack of singly-resonant diagrams. Ultimately, 
the only non-negligible effects are around the peak region, where they are larger than 
5%. The four-lepton transverse momentum is marginally softer with Pythia6 than with 
HERWIG, an effect that is more pronounced in the gg contribution (because of a larger 
amount of radiation there w.r.t. the case of the qq and qg channels, radiation being on 
average softer in Pythia6 than in HERWIG). On the other hand, this larger amount 
of radiation implies that the gg channel is significantly harder than aMC@NLO, and is 
about 10% of the aMC@NLO+(/g result at p T = 0(100 GeV). The scale uncertainty 
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Figure 6: Same as in fig. 4, for the observables introduced in fig. 3. 



shows a behaviour typical of an LO computation, and is vastly different from that of fig. 1; 
it is constant to a good approximation, which is due to a decreasing fj, R dependence with 
increasing p T , compensated by an increasing \x F dependence. 

As was already pointed out in the case of fig. 2, the observables presented in fig. 5 
are quite insensitive to extra radiation, and especially so for the single-inclusive lepton p T . 
Therefore, the similarity between the predictions obtained with HERWIG and Pythia 
has to be expected. To a very good extent, the gg contribution has the same shape as 
the aMC@NLO one. Things are marginally different in the case of the pair transverse 
momentum presented in the right panel; however, any difference between HERWIG and 
Pythia is not significant, being much smaller than the scale uncertainty of the aMC@NLO 
contribution alone (shown in fig. 2). Likewise, the shape change induced in the aMC@NLO 
prediction by the gg contribution is very minor. 

The small-p T dominance for the observables displayed in fig. 6 implies that again the 
HERWIG and Pythia results are very similar. As far as the gg contribution is concerned, 
its shape is almost identical to that of aMC@NLO in the case of Acp. On the other hand, 
the lepton-pair pseudorapidity turns out to be significantly more central than that of the 
qq + qg channels; the impact on the overall aMC@NLO+gg shape is however modest for 
the invariant-mass cuts considered here. 

4. Conclusions 

In this paper, we have used the aMC@NLO framework to study the hadroproduction of 
four leptons, by considering, in particular, the process mediated by the exchange of two 
Z or 7 vector bosons. Di-bosons final states are relevant in both Higgs searches and new- 
physics studies, where precise predictions are an important ingredient to cope with the 
overwhelming irreducible Standard Model backgrounds. We have included in our calcula- 
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tion all contributions which are expected to be important for accurate phenomenological 
analyses, such as singly-resonant diagrams, exact spin correlations and off-shell effects, 
Z/7* interference, and ^-initiated subprocesses. Our computation being fully automated, 
any four-lepton final states can be studied in the very same way, including those that 
feature W/Z interference effects such as e + e~v e v e . 

In any realistic study, precise calculations are as important as the determination of the 
associated theoretical uncertainties. For this reason, we have also given in this work practi- 
cal applications of the fully automatic procedures implemented to this aim in aMC@NLO, 
namely: 

i) The matching of the hard process to both Pythia6 and HERWIG, paving the way 
to systematic studies of the influence of different showers on any observable with 
MC@NLO s ; 

ii) A framework for studying both scale and PDF uncertainties by a simple (and CPU- 
wise costless) reweighting procedure of the unweighted-event samples. 

As for item i), this paper provides the first working example of matching with Pythia6 
for a kinematically non-trivial process. As for item ii), we stress that our procedure allows 
one to obtain precise results with very limited statistics. In fact, the central results and all 
its variations are correlated, since the same events (up to the weights) are used. This is 
manifest in the smoothness of the scale and PDF fractional uncertainty plots presented in 
this paper, which may be compared with, for example, those presented in ref. [18], which 
had been obtained by rerunning the code. 

As far as our phenomenological results are concerned, they can be summarized as 
follows. The impact of NLO corrections on total rates is +40%; for differential distributions, 
they generally correspond to an overall rescaling, but there are cases where they give non- 
trivial kinematics effects. The uncertainties due to PDFs are of the order of 2%. The scale 
dependence we obtain by varying fj, R and fi F independently is also of the order of 2% for 
the qq/qg channels, and of the order of 20% in the case of the ^-initiated subprocess. 
Parton-level NLO and fully showered NLO distributions are generally in good agreement, 
except in the very few cases in which the former is not expected to give sensible predictions. 
These findings are rather independent of possible Z-identification cuts applied to enhance 
doubly-resonant contributions. Generally speaking, tiny differences are observed between 
HERWIG and Pythia, with the exception of distributions dominated by resummation 
effects, where they are more pronounced. The contribution of the gg channel is of the order 
of 5% of the total with a cut of 30 GeV applied to unlike-sign lepton pairs. 

We have explicitly verified that the inclusion of the gg —> H — >■ ZZ contribution 
(not considered for producing the results presented here) is straightforward, and does not 
increase significantly the running time. This opens the way to simulating Higgs production 
with four-lepton decay channels by including the signal/background interference. Clearly, 
for this to be phenomenologically sensible, NLO corrections to the gg — > H signal should 

3 The matching with Herwig++ and Pythia8 will become available as well. 
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to be included as well, which could be achieved in several ways, with different levels of 
approximation and automation [27, 28, 29]. 

As a last remark, we point out that ready-to-shower four-lepton NLO (i.e., result- 
ing from qq /qg-imtiaied processes) event files are available at the aMC@NLO web page 
http://amcatnlo.cern.ch, for all possible (massless) lepton flavour combinations 4 . 
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A. Definition of W coefficients for NLO cross sections 

We collect in this section the expressions of the coefficients that enter eq. (2.17). These 
are taken from ref. [23], whose equation (n) will be denoted by eq. (I.n) here; the reader 
can find fuller details in that paper. The master equation for the determination of all the 
coefficients is eq. (I.C.14). 



• Event 

Only the (n + l)-body cross section contributes to the event. Therefore, from eq. (1.4.29) 
and eq. (1.6.7) one gets: 

W { E) = w real (/C£\) , (A.l) 

W { F E) = , (A.2) 
t?f } = 0, (A.3) 



where 



% Sjj(/C n +l) T f\ A\ 

WrealK^n+l) ~ — 2 h+2 / T J Bj > l A - 4 J 



6(1 -wy) g? +2 (nn) 



S« = ((1 - 1*)&M (B+1 ' 0) ) Sy^* (B+1) - (A.5) 



M 

with |>( n+1 ) = $( n+1 )/&. The quantity defined in eq. (A.5) is the combination of the real- 
emission matrix element, prefactors, and the phase-space factor that is routinely used in 
the FKS subtraction. We point out that according to our conventions the matrix elements 



4 Processes that feature two W propagators interfere at the NLO with tW production. We have excluded 
the latter contribution by setting the 6-quark PDF equal to zero. 
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include the coupling constant, hence eq. (A. 4) is independent of g s , since we understand 
that the factorization scale \i R there is computed according to the kinematic configuration 
/C n+ i. The coupling-constant dependence has been factored out in eq. (2.17). 



• Collinear counterevent 

The coefficient receives contributions from the (n + l)-body and the degenerate 

(n + l)-body cross sections. The latter can be read from eq. (1.4.41) (to be definite, we 
assume that the FKS sister is the parton coming from the left). One obtains: 



W, 



(CO 



(CO 



w. 



+ 



a s {^R ] ) J -5(0) 

i©t 




i" log ^ +2 — 



(S) 



1 



-^• 0) (4 c + \) J -j^w^ (4?r) ^ , 



(A.6) 



(A.7) 
(A.8) 



The normalization factor in eqs. (A.6) and (A.7) is due to the fact that the degenerate 
(n + l)-body contribution is integrated together with the pure (n + l)-body one, but is 
originally defined with a quasi- n-body measure (see eq. (1.4.41)). 



• Soft counterevent, n-body contribution, and Born 

The coefficient receives contributions from the (n + l)-body and the n-body cross sec- 
tions. The latter is the sum of the Born (eq. (1.4.4)), collinear (eq. (1.4.5)), soft (eq. (1.4.12)) 
and finite virtual (eq. (1.4.14)) contributions, plus the last term on the r.h.s. of eq. (I.C.14), 
which we shall call the RGE term here. The expression for the finite part of the virtual 
contribution that includes the full scale dependence must be read from eq. (I.C.13). In 
ref. [23] we have used the Ellis-Sexton scale wherever possible, which now renders it easy to 
extract the factorization and renormalization scale dependences - more specifically, these 
dependences are present only in one term in eq. (1.4.5), in one term in eq. (I.C.13), and in 
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the RGE term. This implies: 



W, 



(S) _ 



(4 s ?,) 



bt {S)\ 

g s {fl R ) („,!) / 2 2 n 2. r (S) \ 



+ 



1 



26+2 /„ (5) x 27T 



:■(<?) 



(n,0) / r (S) \ 

kl y^n+i J 



(/Cgj J Bi , (A.9) 



x ( 7 (X!) + 2C(Xi) log£ cut + 7 (X 2 ) + 2C(X 2 ) log£ CMi ) , 



(A.10) 



1 



«s(/4f } ) 



2tt/3 & 



(S) 



xM™ -^^ 



(B) 



(A.11) 
(A.12) 



The normalization factor is the analogue of that introduced in eq. (A. 6), and its form 
can be found e.g. in eq. (1.6.14) (note that it includes an Sij factor). In eq. (A.9), by- 
setting fi = Q one eliminates the first term on the r.h.s. of eq. (1.4.6); this is included here 
in eq. (A.10). 

• Soft-collinear counterevent 

Finally, the results for W^ sc ^ can be obtained by computing the soft limit of W^ c \ and by 
taking into account the fact that the original equations contained plus distributions. One 
obtains: 
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